function obj=exp_u(y, rho, kappa, zeta, logit_coef, X)


X_new=[y'; X]; 

temp=prod(y.^zeta,2);

prob=exp(sum(logit_coef.*X_new,1))/(1+exp(sum(logit_coef.*X_new,1)));

obj=prob*u(rho,temp-kappa)+(1-prob)*u(rho, temp);


end